Multi-omics analysis for samples with matched brain and retina tissue from same subjects¶

  • GLDS-203 v6 and GLDS-203 v8
  • Select subjects that have methylation/RNA-Seq from matched tissue
  • Define promoter regions as +/- 2kb from TSS
  • Select all CpG sites within promoter of each gene
  • Define gene methylation level as median methylation within promoter of each gene for each sample
  • Option 1: For each group, calculate correlation between RNA-Seq and RRBS counts per gene for each tissue separarely and then intersect between tissues
  • Option 2: For each group, calculate correlation between RNA-Seq or RRBS counts per gene for each assay separately and then intersect between assays
In [145]:
#Load libraries

#General
library(deconvSeq)
library(data.table)
library(plyr)
library(dplyr)
library(tidyverse)
library(ggplot2)
library(DESeq2)
library(matrixStats)
library(ComplexUpset)
library("RColorBrewer")
library(treemapify)
library(ggtree)
library(cowplot)
library(pheatmap)
library(MDSeq)
library(matrixStats)

#Annotation
organism="org.Mm.eg.db"
library(organism, character.only = TRUE)
library("genomation")
library("GenomicRanges")
library("biomaRt")

#GSEA and ORA
library(clusterProfiler)
library(enrichplot)
library(ggnewscale)
library(rrvgo)
In [146]:
#Get gene mappings and biotypes
httr::set_config(httr::config(ssl_verifypeer = FALSE))
mart <- useMart('ENSEMBL_MART_ENSEMBL')
mart <- useDataset('mmusculus_gene_ensembl', mart)
annotLookup <- getBM(
  mart = mart,
  attributes = c(
    'ensembl_gene_id',
    'ensembl_transcript_id',
    'external_gene_name',
    'gene_biotype',
    'entrezgene_id'),
  uniqueRows = TRUE)
In [ ]:
#Read in transcript and CpG island annotations
gene.obj=readTranscriptFeatures("/home/jupyter/MOUNT/references/mm10.GRCm38.96.bed", up.flank = 2000,down.flank = 2000, remove.unusual=TRUE, unique.prom = TRUE)
cpg.obj=readFeatureFlank("/home/jupyter/MOUNT/references/cpgi.mm10.bed",feature.flank.name=c("CpGi","shores"), remove.unusual = TRUE,)
bedfile<-fread("/home/jupyter/MOUNT/references/mm10.GRCm38.96.bed",header=F)
bedfile<-subset(bedfile,select=c("V1","V4","V6","V2","V3"))
colnames(bedfile)<-c("chr","feature.name","transStrand","transStart","transStop")
In [172]:
##Define functions
#Function for annotating methylation matrix
annotateMeth <- function(data){

    gr <- makeGRangesFromDataFrame(data, keep.extra.columns=FALSE,
                                      start.field="start",
                                      end.field="end")

    cpgAnn=annotateWithFeatureFlank(gr,
                                    cpg.obj$CpGi,cpg.obj$shores,feature.name="CpGi",flank.name="shores")
    cpgMat<-as.data.frame(getMembers(cpgAnn))
    cpgMat$CpG<-ifelse(cpgMat$CpGi==1,"Island",ifelse(cpgMat$shores==1,"Shore","Other"))
    data<-cbind(data,cpgMat)
        
    #Gene annotation
    geneAnn=annotateWithGeneParts(gr,gene.obj)
    geneMat<-cbind(getAssociationWithTSS(geneAnn), as.data.frame(getMembers(geneAnn)))
    data<-cbind(data,geneMat)
    data<-merge(data,bedfile,by=c("chr","feature.name"))
    #Workaround for a bug in genomation (https://github.com/BIMSBbioinfo/genomation/issues/203)
    data$intron<-ifelse((data$start<=data$transStart & data$transStrand=="+"),0,
                    ifelse((data$start>=data$transStop & data$transStrand=="-"),0,data$intron))
    data$Region<-ifelse(data$prom==1,"Promoter",ifelse(data$exon==1 | data$intron==1,"Gene body","Intergenic"))
    #Filter on genomic context and significance
    data<-subset(data,Region=="Promoter")
    data<-merge(data,annotLookup,by.x=c("feature.name"),by.y=c("ensembl_transcript_id"))
    data
}

#Function for upset plot
createUpsetPlot <- function(data) {
    df<-reshape2::dcast(data, gene~Type, length)
    Groups = colnames(df)[2:3]
    plotUpset<-upset(df, Groups, name='Groups', width_ratio=0.1, min_size=0, min_degree=1, set_sizes=FALSE,
      queries=list(
        upset_query(set='Control', fill='gray'),
        upset_query(set='IR', fill='magenta')
      ),
      themes=upset_default_themes(text=element_text(size=24,face="bold"), 
                                  axis.title=element_text(face="bold",size=26),
                                  panel.border = element_rect(colour = "black", fill=NA, size=1)),
      sort_intersections_by='cardinality',
      base_annotations=list('Size'=(intersection_size(counts=TRUE))),
       matrix=intersection_matrix(geom=geom_point(shape='square filled',size=5)))
     plotUpset
}
In [141]:
#Sample manifest
manifest<-fread("/home/jupyter/MOUNT/integrated/manifest.txt",header=T,sep="\t")
manifest<- manifest %>% dplyr::filter(Brain_rna!="Mmus_C57_6J_BRN_HLU_IRC_4mon_Rep5_M93") %>% dplyr::filter(Tissue=="Both" & Group %in% c("IR","Control"))
manifest[1:2,]
A data.table: 2 × 13
SubjectRepTissueSampleTimepointGroupBrain_methRetina_methBrain_rnaRetina_rnaCountBrain_ageRetina_age
<chr><chr><chr><chr><chr><chr><chr><chr><chr><chr><int><dbl><dbl>
M83Rep1BothHLLC_IRC_4mon_Rep1_M834mControlCFG1778CFG1837Mmus_C57_6J_BRN_HLLC_IRC_4mon_Rep1_M83Mmus_C57_6J_RTN_HLLC_IRC_4mon_Rep1_M83413.33.4
M85Rep2BothHLLC_IRC_4mon_Rep2_M854mControlCFG1780CFG1839Mmus_C57_6J_BRN_HLLC_IRC_4mon_Rep2_M85Mmus_C57_6J_RTN_HLLC_IRC_4mon_Rep2_M85414.1 NA
In [ ]:
#Process and merge brain methylation data
samples_202=manifest$Brain_meth
files_202=paste("GLDS-202_Epigenomics_R1_",samples_202,
                "_trimmed.fq_trimmed_bismark_bt2.bismark.cov.gz",sep="")
setwd("/home/jupyter/MOUNT/integrated/bismark/glds202")
methmat_202_all = getmethmat(filnames=files_202, sample.id=samples_202, filtype="bismarkCoverage")
methmat_202_all<-annotateMeth(methmat_202_all)
In [ ]:
#Process and merge retina methylation data
samples_203=manifest$Retina_meth
files_203=paste("GLDS-203_Epigenomics_R1_",samples_203,
                "_trimmed.fq_trimmed_bismark_bt2.bismark.cov.gz",sep="")
setwd("/home/jupyter/MOUNT/integrated/bismark/glds203")
methmat_203_all = getmethmat(filnames=files_203, sample.id=samples_203, filtype="bismarkCoverage")
methmat_203_all<-annotateMeth(methmat_203_all)
In [150]:
manifest
A data.table: 9 × 13
SubjectRepTissueSampleTimepointGroupBrain_methRetina_methBrain_rnaRetina_rnaCountBrain_ageRetina_age
<chr><chr><chr><chr><chr><chr><chr><chr><chr><chr><int><dbl><dbl>
M83Rep1BothHLLC_IRC_4mon_Rep1_M834mControlCFG1778CFG1837Mmus_C57_6J_BRN_HLLC_IRC_4mon_Rep1_M83Mmus_C57_6J_RTN_HLLC_IRC_4mon_Rep1_M83413.33.4
M85Rep2BothHLLC_IRC_4mon_Rep2_M854mControlCFG1780CFG1839Mmus_C57_6J_BRN_HLLC_IRC_4mon_Rep2_M85Mmus_C57_6J_RTN_HLLC_IRC_4mon_Rep2_M85414.1 NA
M75Rep1BothHLLC_IR_4mon_Rep1_M75 4mIR CFG1771CFG1831Mmus_C57_6J_BRN_HLLC_IR_4mon_Rep1_M75 Mmus_C57_6J_RTN_HLLC_IR_4mon_Rep1_M75 415.32.3
M80Rep2BothHLLC_IR_4mon_Rep2_M80 4mIR CFG1775CFG1835Mmus_C57_6J_BRN_HLLC_IR_4mon_Rep2_M80 Mmus_C57_6J_RTN_HLLC_IR_4mon_Rep2_M80 4 NA3.7
M82Rep4BothHLLC_IR_4mon_Rep4_M82 4mIR CFG1777CFG1836Mmus_C57_6J_BRN_HLLC_IR_4mon_Rep4_M82 Mmus_C57_6J_RTN_HLLC_IR_4mon_Rep3_M82 414.31.4
M84Rep5BothHLLC_IR_4mon_Rep5_M84 4mIR CFG1779CFG1838Mmus_C57_6J_BRN_HLLC_IR_4mon_Rep5_M84 Mmus_C57_6J_RTN_HLLC_IR_4mon_Rep4_M84 4 NA2.4
M89Rep6BothHLLC_IR_4mon_Rep6_M89 4mIR CFG1783CFG1842Mmus_C57_6J_BRN_HLLC_IR_4mon_Rep6_M89 Mmus_C57_6J_RTN_HLLC_IR_4mon_Rep5_M89 411.92.8
M95Rep3BothHLLC_IRC_4mon_Rep3_M954mControlCFG1788CFG1846Mmus_C57_6J_BRN_HLLC_IRC_4mon_Rep3_M95Mmus_C57_6J_RTN_HLLC_IRC_4mon_Rep3_M95412.82.2
M97Rep4BothHLLC_IRC_4mon_Rep4_M974mControlCFG1789CFG1847Mmus_C57_6J_BRN_HLLC_IRC_4mon_Rep4_M97Mmus_C57_6J_RTN_HLLC_IRC_4mon_Rep4_M97414.22.5
In [151]:
#Only keep sites with >=10 reads in all samples
selectedCols <- methmat_202_all[grep("coverage", names(methmat_202_all), fixed = T)]
methmat_202_all <- methmat_202_all[apply(selectedCols, 1, 
                                   function(x) all(!is.na(x) & x>=10)), ]
selectedCols <- methmat_203_all[grep("coverage", names(methmat_203_all), fixed = T)]
methmat_203_all <- methmat_203_all[apply(selectedCols, 1, 
                                   function(x) all(!is.na(x) & x>=10)), ]
colMeans(methmat_202_all[grep("coverage", names(methmat_202_all), fixed = T)])
colMeans(methmat_203_all[grep("coverage", names(methmat_203_all), fixed = T)])
dim(methmat_202_all)
coverage1
37.1496370786918
coverage2
41.7942083394515
coverage3
36.2953096042193
coverage4
58.602116015763
coverage5
35.5482459795568
coverage6
123.85893492815
coverage7
41.5014809200092
coverage8
45.5277340138025
coverage9
24.235589967048
coverage1
35.9773426829348
coverage2
74.5546402734301
coverage3
20.0632029355832
coverage4
36.3882132284226
coverage5
35.5660375508924
coverage6
66.8226217005731
coverage7
40.9288225645181
coverage8
35.0406436455604
coverage9
43.0981250694358
  1. 755949
  2. 49
In [152]:
#Calculate methylation % for all sites/samples
methmat_202<-methmat_202_all
methmat_203<-methmat_203_all
rm(methmat_202_all)
rm(methmat_203_all)

for (i in 1:9){
    numCs <- rlang::sym(paste("numCs", i, sep = ""))
    coverage <- rlang::sym(paste("coverage", i, sep = ""))
    methmat_202<-methmat_202 %>% dplyr::mutate(!!paste0("ratio_",manifest$Brain_meth[i]) := !!numCs/!!coverage)
    methmat_203<-methmat_203 %>% dplyr::mutate(!!paste0("ratio_",manifest$Retina_meth[i]) := !!numCs/!!coverage)
}

#Calculate ratios (per site) and filter sites with sd<0.02
methmat_202 <- methmat_202 %>% mutate(mean = select(., starts_with("ratio_")) %>% rowMeans(na.rm = TRUE),
                                     sd = rowSds(as.matrix(select(.,starts_with("ratio_"))),na.rm = TRUE)) %>%
                               filter(sd > 0.02)
methmat_203 <- methmat_203 %>% mutate(mean = select(., starts_with("ratio_")) %>% rowMeans(na.rm = TRUE),
                                     sd = rowSds(as.matrix(select(.,starts_with("ratio_"))),na.rm = TRUE)) %>%
                               filter(sd > 0.02)
dim(methmat_202)
dim(methmat_203)
  1. 187432
  2. 60
  1. 168854
  2. 60
In [153]:
for (i in 1:9){
    numCs <- rlang::sym(paste("numCs", i, sep = ""))
    coverage <- rlang::sym(paste("coverage", i, sep = ""))
    methmat_202<-methmat_202 %>% dplyr::mutate(!!coverage := ifelse(!!sym(numCs)==0,NA,!!sym(coverage)), !!numCs := ifelse(!!sym(numCs)==0,NA,!!sym(numCs)))
    methmat_203<-methmat_203 %>% dplyr::mutate(!!coverage := ifelse(!!sym(numCs)==0,NA,!!sym(coverage)), !!numCs := ifelse(!!sym(numCs)==0,NA,!!sym(numCs)))
}
In [154]:
#Calculate weighted methylation for each gene
methmat_202_counts<- methmat_202 %>% dplyr::select(c("external_gene_name",contains(c("numCs","coverage")))) %>%
                    dplyr::group_by(external_gene_name) %>%  
                    summarise(across(starts_with(c('numCs', 'coverage')), sum, na.rm=TRUE))
methmat_203_counts<- methmat_203 %>% dplyr::select(c("external_gene_name",contains(c("numCs","coverage")))) %>%
                    dplyr::group_by(external_gene_name) %>%  
                    summarise(across(starts_with(c('numCs', 'coverage')), sum, na.rm=TRUE))
for (i in 1:9){
    coverage <- rlang::sym(paste("coverage", i, sep = ""))
    numCs <- rlang::sym(paste("numCs", i, sep = ""))
    methmat_202_counts<-methmat_202_counts %>% mutate(!!paste("ratio",manifest$Brain_meth[i],sep="_") := !!numCs/!!coverage)
    methmat_203_counts<-methmat_203_counts %>% mutate(!!paste("ratio",manifest$Retina_meth[i],sep="_") := !!numCs/!!coverage)
}

rm(methmat_202)
rm(methmat_203)
In [155]:
#Methylation ratio distribution
df_202 <- as.data.frame(gather(as.data.frame(methmat_202_counts %>% select(contains("ratio_")))))
df_202$key<-gsub("ratio_","",df_202$key)
df_203 <- as.data.frame(gather(as.data.frame(methmat_203_counts %>% select(contains("ratio_")))))
df_203$key<-gsub("ratio_","",df_203$key)

df_202<-merge(df_202,manifest,by.x=c("key"),by.y=c("Brain_meth"))
df_203<-merge(df_203,manifest,by.x=c("key"),by.y=c("Retina_meth"))

ggplot(df_202) + geom_density(aes(x=value,color=Sample))
ggplot(df_203) + geom_density(aes(x=value,color=Sample))

Warning message:
“Removed 21351 rows containing non-finite values (stat_density).”
Warning message:
“Removed 19092 rows containing non-finite values (stat_density).”
In [158]:
#Compute PCs
tmp<-methmat_203_counts[,grep("ratio", names(methmat_203_counts))]
tmp[is.na(tmp)] <- 0
dfWide.pca <- prcomp(t(tmp), center = TRUE, scale. = TRUE)
dfWide.pca <- as.data.frame(dfWide.pca$x)
In [162]:
#PCA plots
#dfWide.pca$Tissue=rep(c("BRN","RTN"),each=length(manifest$Sample))
dfWide.pca$Group=rep(manifest$Group,times=1)
ggplot(
  dfWide.pca,
  aes(x = PC3, y = PC4,color=Group)
  ) +
  scale_shape_manual(values=c(17,16)) +
  geom_point(size=6, alpha = 0.7)+
  geom_hline(yintercept = 0, colour = "gray65") +
  geom_vline(xintercept = 0, colour = "gray65") + theme_classic()

Process RNA-seq data¶

In [163]:
genesCoding<-unique(as.data.frame(subset(annotLookup, gene_biotype=="protein_coding" | grepl("RNA",gene_biotype), select=c("ensembl_gene_id"))))
colnames(genesCoding) <-c("gene")

#BRN dataset (GLDS-202)
counts202<-read.delim("/home/jupyter/MOUNT/integrated/GLDS-202_rna_seq_Unnormalized_Counts.csv",sep=",",row.names=1,header=TRUE)
counts202<-counts202 %>% dplyr::filter(!grepl("ERCC",row.names(counts202))) %>% 
            dplyr::select(manifest$Brain_rna)
genesCounts<-as.data.frame(row.names(counts202))
colnames(genesCounts) <- c("gene")
merged<-merge(genesCounts,genesCoding,by=c("gene"))
counts202 <- counts202[merged$gene,]

coldata202<-as.data.frame(colnames(counts202))
colnames(coldata202)<-c("SampleName")
coldata202$Group<-manifest$Group
dds202 <- DESeqDataSetFromMatrix(countData = round(counts202),colData = coldata202,design =~ 1 )
ddsNofilt202 <- estimateSizeFactors(dds202)
keep <- rowSums(counts(ddsNofilt202, normalized=TRUE) >= 5 ) >= 4 & rowSums(counts(ddsNofilt202, normalized=TRUE) == 0 ) <= 4
dds202 <- ddsNofilt202[keep,]
counts202<-log2(as.data.frame(counts(dds202, normalized=TRUE))+1)
counts202$ensembl_gene_id<-rownames(counts202)
counts202<-merge(counts202,unique(subset(annotLookup,select=c("ensembl_gene_id","external_gene_name"))),by=c("ensembl_gene_id"))

#RTN dataset (GLDS-203)
counts203<-read.delim("/home/jupyter/MOUNT/integrated/GLDS-203_rna_seq_Unnormalized_Counts.csv",sep=",",row.names=1,header=TRUE)
counts203<-counts203 %>% filter(!grepl("ERCC",row.names(counts203))) %>%
            dplyr::select(manifest$Retina_rna) 
genesCounts<-as.data.frame(row.names(counts203))
colnames(genesCounts) <-c("gene")
merged<-merge(genesCounts,genesCoding,by=c("gene"))
counts203 <- counts203[merged$gene,]
coldata203<-as.data.frame(colnames(counts203))
colnames(coldata203)<-c("SampleName")
coldata203$Group<-manifest$Group
dds203 <- DESeqDataSetFromMatrix(countData = round(counts203),colData = coldata203,design =~ 1 )

ddsNofilt203 <- estimateSizeFactors(dds203)
keep <- rowSums(counts(ddsNofilt203, normalized=TRUE) >= 5 ) >= 4 & rowSums(counts(ddsNofilt203, normalized=TRUE) == 0 ) <= 4
dds203 <- ddsNofilt203[keep,]
counts203<-log2(as.data.frame(counts(dds203, normalized=TRUE))+1)
counts203$ensembl_gene_id<-rownames(counts203)
counts203<-merge(counts203,unique(subset(annotLookup,select=c("ensembl_gene_id","external_gene_name"))),by=c("ensembl_gene_id"))
dim(counts202)
dim(counts203)

converting counts to integer mode

converting counts to integer mode
  1. 15941
  2. 11
  1. 16179
  2. 11

Within assay: Correlation between brain and retina in each group¶

Within each exposure group, calculate Pearson correlation for each gene promoter (CpGi/s) using brain-retina matched pairs

RRBS¶

In [165]:
#Merge brain and retina data on common genes
mergedMeth<-merge(methmat_202_counts %>% dplyr::select(c("external_gene_name",contains("ratio"))),
                  methmat_203_counts %>% dplyr::select(c("external_gene_name",contains("ratio"))),
                  by=c("external_gene_name"))
mergedMeth[is.na(mergedMeth)] <- 0
In [168]:
#Calculate correlation between brain and retina for each gene within each group
for (group in c("IR","Control")){
    tmp202<-as.matrix(mergedMeth[,paste0("ratio_",subset(manifest,Group==group)$Brain_meth)])
    tmp203<-as.matrix(mergedMeth[,paste0("ratio_",subset(manifest,Group==group)$Retina_meth)])
    corMatTmp<-suppressWarnings(sapply(seq.int(dim(tmp202)[1]), 
           function(i) cor.test(tmp202[i,], tmp203[i,])))
    corMatTmp<-as.data.frame(t(corMatTmp)) %>% dplyr::select(statistic,p.value,estimate,conf.int) %>% 
                                    rename_with( ~ paste0(.x,".",group))
    assign(paste0("corMat_",group),corMatTmp)
}
corMatMeth<-cbind(corMat_IR,corMat_Control)
corMatMeth$gene<-mergedMeth$external_gene_name
In [171]:
posCorMeth[1,]
A data.frame: 1 × 2
geneType
<chr><chr>
11110006O24RikIR
In [173]:
#Positively correlated
posCorMeth<-unique(rbind(as.data.frame(corMatMeth %>% dplyr::filter(p.value.IR<0.05 & estimate.IR>0) 
                            %>% dplyr::select(gene) %>% mutate(Type="IR")),
              as.data.frame(corMatMeth %>% dplyr::filter(p.value.Control<0.05 & estimate.Control>0) 
                            %>% dplyr::select(gene) %>% mutate(Type="Control"))))

#Negatively correlated
negCorMeth<-unique(rbind(as.data.frame(corMatMeth %>% dplyr::filter(p.value.IR<0.05 & estimate.IR<0) 
                            %>% dplyr::select(gene) %>% mutate(Type="IR")),
              as.data.frame(corMatMeth %>% dplyr::filter(p.value.Control<0.05 & estimate.Control<0) 
                            %>% dplyr::select(gene) %>% mutate(Type="Control"))))
options(repr.plot.width=6, repr.plot.height=8)
plot_grid(createUpsetPlot(posCorMeth), createUpsetPlot(negCorMeth),ncol=1)
posCorMeth$assay<-"RRBS"
negCorMeth$assay<-"RRBS"

Using Type as value column: use value.var to override.

Using Type as value column: use value.var to override.

RNA-Seq¶

In [174]:
mergedRna<-merge(counts202,counts203,by="ensembl_gene_id")
for (group in c("IR","Control")){
    tmp202<-as.matrix(t(apply(mergedRna[,subset(manifest,Group==group)$Brain_rna],1,scale)))
    tmp203<-as.matrix(t(apply(mergedRna[,subset(manifest,Group==group)$Retina_rna],1,scale)))
    corMatTmp<-suppressWarnings(sapply(seq.int(dim(tmp202)[1]), 
           function(i) cor.test(tmp202[i,], tmp203[i,])))
    corMatTmp<-as.data.frame(t(corMatTmp)) %>% dplyr::select(statistic,p.value,estimate,conf.int) %>% 
                                    rename_with( ~ paste0(.x,".",group))
    assign(paste0("corMatRna_",group),corMatTmp)
}

corMatRna<-cbind(corMatRna_IR,corMatRna_Control)
corMatRna$ensembl_gene_id<-mergedRna$ensembl_gene_id
corMatRna<-merge(corMatRna,unique(subset(annotLookup,select=c("ensembl_gene_id","external_gene_name"))),by=c("ensembl_gene_id"))
corMatRna <- corMatRna %>% rename(gene=external_gene_name)
In [175]:
#Positively correlated
posCorRna<-unique(rbind(as.data.frame(corMatRna %>% dplyr::filter(p.value.IR<0.05 & estimate.IR>0) 
                            %>% dplyr::select(gene) %>% mutate(Type="IR")),
              as.data.frame(corMatRna %>% dplyr::filter(p.value.Control<0.05 & estimate.Control>0) 
                            %>% dplyr::select(gene) %>% mutate(Type="Control"))))

#Negatively correlated
negCorRna<-unique(rbind(as.data.frame(corMatRna %>% dplyr::filter(p.value.IR<0.05 & estimate.IR<0) 
                            %>% dplyr::select(gene) %>% mutate(Type="IR")),
              as.data.frame(corMatRna %>% dplyr::filter(p.value.Control<0.05 & estimate.Control<0) 
                            %>% dplyr::select(gene) %>% mutate(Type="Control"))))
options(repr.plot.width=6, repr.plot.height=8)
plot_grid(createUpsetPlot(posCorRna), createUpsetPlot(negCorRna),nrow=2)
posCorRna$assay<-"RNA"
negCorRna$assay<-"RNA"

Using Type as value column: use value.var to override.

Using Type as value column: use value.var to override.
In [181]:
allPosCor<-rbind(posCorMeth,posCorRna)
allPosCor$Type<-paste(allPosCor$assay,allPosCor$Type,sep="_")

    df<-reshape2::dcast(allPosCor, gene~Type, length)
    df[1,]
    Groups = colnames(df)[2:5]
    plotUpset<-upset(df, Groups, name='Groups', width_ratio=0.1, min_size=0, min_degree=1, set_sizes=FALSE,
      queries=list(
        upset_query(set='RNA_Control', fill='lightblue'),
        upset_query(set='RNA_IR', fill='orange'),
        upset_query(set='RRBS_Control', fill='blue'),
        upset_query(set='RRBS_IR', fill='brown')
      ),
      themes=upset_default_themes(text=element_text(size=24,face="bold"), 
                                  axis.title=element_text(face="bold",size=26),
                                  panel.border = element_rect(colour = "black", fill=NA, size=1)),
      sort_intersections_by='cardinality',
      base_annotations=list('Size'=(intersection_size(counts=TRUE))),
       matrix=intersection_matrix(geom=geom_point(shape='square filled',size=5)))
options(repr.plot.width=15, repr.plot.height=7)
     plotUpset
dfSum = df %>% mutate(sum = select(., !starts_with("gene") & !contains("Control")) %>% rowSums(na.rm = TRUE))

Using assay as value column: use value.var to override.
A data.frame: 1 × 5
geneRNA_ControlRNA_IRRRBS_ControlRRBS_IR
<chr><int><int><int><int>
10610010K14Rik0010
In [189]:
subset(dfSum, sum>0 & RRBS_Control==0 & RNA_Control==0 & RNA_IR==1 & RRBS_IR==1)$gene
  1. 'Agbl2'
  2. 'Begain'
  3. 'Dlg4'
  4. 'Elfn1'
  5. 'Eps8l1'
  6. 'Fam89b'
  7. 'Fbxl7'
  8. 'Gm28778'
  9. 'Gpr89'
  10. 'Grm8'
  11. 'Kalrn'
  12. 'Ldlrad3'
  13. 'Mmp15'
  14. 'Nog'
  15. 'Pdgfra'
  16. 'Pkdcc'
  17. 'Rasgrp4'
  18. 'Rpl8'
  19. 'Slc18b1'
  20. 'Taok1'
  21. 'Vipr1'
  22. 'Ybx1'

Within tissue: Correlation between RNA-seq and RRBS in each group¶

Within each exposure group for each tissue type, calculate Pearson correlation for each gene promoter (CpGi/s) using RRBS-RNA-seq matched pairs

In [190]:
methmat_203_counts[is.na(methmat_203_counts)] <- 0
methmat_202_counts[is.na(methmat_202_counts)] <- 0
merged203<-merge(counts203,methmat_203_counts,by="external_gene_name")
for (group in c("IR","Control")){
    tmpMeth<-as.matrix(merged203[,paste0("ratio_",subset(manifest,Group==group)$Retina_meth)])
    tmpRna<-as.matrix(t(apply(merged203[,subset(manifest,Group==group)$Retina_rna],1,scale)))
    corMat<-suppressWarnings(sapply(seq.int(dim(tmpMeth)[1]), 
           function(i) cor.test(tmpMeth[i,], tmpRna[i,])))
    corMat<-as.data.frame(t(corMat)) %>% dplyr::select(statistic,p.value,estimate,conf.int) %>% 
                                    rename_with( ~ paste0(.x,".",group))
    assign(paste0("corMatRetina_",group),corMat)
}
                                    
merged202<-merge(counts202,methmat_202_counts,by="external_gene_name")
for (group in c("IR","Control")){
    samples_rna<-subset(manifest,Group==group)$Brain_rna
    samples_meth<-subset(manifest,Group==group)$Brain_meth
    tmpMeth<-as.matrix(merged202[,paste0("ratio_",samples_meth)])
    tmpRna<-as.matrix(t(apply(merged202[,samples_rna],1,scale)))
    corMat<-suppressWarnings(sapply(seq.int(dim(tmpMeth)[1]), 
           function(i) cor.test(tmpMeth[i,], tmpRna[i,])))
    corMat<-as.data.frame(t(corMat)) %>% dplyr::select(statistic,p.value,estimate,conf.int) %>% 
                                    rename_with( ~ paste0(.x,".",group))
    assign(paste0("corMatBrain_",group),corMat)
}
In [246]:
samples_rna<-subset(manifest,Group==group)$Brain_rna
samples_meth<-subset(manifest,Group==group)$Brain_meth

ggplot(subset(merged202, !!subset(manifest,Group==group)$Brain_rna[2] > 5),aes(x=get(subset(manifest,Group==group)$Brain_rna[2]), 
                                   y=get(paste0("ratio_",samples_meth[2])))) + geom_point() + geom_smooth(method="lm")

`geom_smooth()` using formula 'y ~ x'
In [191]:
corMatBrain<-cbind(corMatBrain_IR,corMatBrain_Control)
corMatRetina<-cbind(corMatRetina_IR,corMatRetina_Control)
corMatBrain$gene<-merged202$external_gene_name
corMatRetina$gene<-merged203$external_gene_name
dim(corMatBrain)
dim(corMatRetina)
  1. 12385
  2. 9
  1. 12085
  2. 9
In [229]:
t(subset(merged202, external_gene_name=="Fktn"))
A matrix: 38 × 1 of type chr
3864
external_gene_nameFktn
ensembl_gene_idENSMUSG00000028414
Mmus_C57_6J_BRN_HLLC_IRC_4mon_Rep1_M839.753355
Mmus_C57_6J_BRN_HLLC_IRC_4mon_Rep2_M859.834692
Mmus_C57_6J_BRN_HLLC_IR_4mon_Rep1_M759.672554
Mmus_C57_6J_BRN_HLLC_IR_4mon_Rep2_M809.708237
Mmus_C57_6J_BRN_HLLC_IR_4mon_Rep4_M829.369445
Mmus_C57_6J_BRN_HLLC_IR_4mon_Rep5_M849.746524
Mmus_C57_6J_BRN_HLLC_IR_4mon_Rep6_M899.762727
Mmus_C57_6J_BRN_HLLC_IRC_4mon_Rep3_M9510.03492
Mmus_C57_6J_BRN_HLLC_IRC_4mon_Rep4_M979.790137
numCs12
numCs23
numCs31
numCs474
numCs567
numCs6118
numCs75
numCs85
numCs91
coverage1106
coverage2115
coverage359
coverage41350
coverage5555
coverage62688
coverage7220
coverage8240
coverage940
ratio_CFG17780.01886792
ratio_CFG17800.02608696
ratio_CFG17710.01694915
ratio_CFG17750.05481481
ratio_CFG17770.1207207
ratio_CFG17790.04389881
ratio_CFG17830.02272727
ratio_CFG17880.02083333
ratio_CFG17890.025
In [192]:
#Positively correlated between RRBS and RNA-Seq
posCorRetina<-unique(rbind(as.data.frame(corMatRetina %>% dplyr::filter(p.value.IR<0.05 & estimate.IR>0) 
                            %>% dplyr::select(gene) %>% mutate(Type="IR")),
              as.data.frame(corMatRetina %>% dplyr::filter(p.value.Control<0.05 & estimate.Control>0) 
                            %>% dplyr::select(gene) %>% mutate(Type="Control"))))

#Negatively correlated between RRBS and RNA-Seq
negCorRetina<-unique(rbind(as.data.frame(corMatRetina %>% dplyr::filter(p.value.IR<0.05 & estimate.IR<0) 
                            %>% dplyr::select(gene) %>% mutate(Type="IR")),
              as.data.frame(corMatRetina %>% dplyr::filter(p.value.Control<0.05 & estimate.Control<0) 
                            %>% dplyr::select(gene) %>% mutate(Type="Control"))))
options(repr.plot.width=6, repr.plot.height=8)
plot_grid(createUpsetPlot(posCorRetina), createUpsetPlot(negCorRetina),ncol=1)
posCorRetina$tissue<-"RTN"
negCorRetina$tissue<-"RTN"

Using Type as value column: use value.var to override.

Using Type as value column: use value.var to override.
In [193]:
#Positively correlated between RRBS and RNA-Seq
posCorBrain<-unique(rbind(as.data.frame(corMatBrain %>% dplyr::filter(p.value.IR<0.05 & estimate.IR>0) 
                            %>% dplyr::select(gene) %>% mutate(Type="IR")),
              as.data.frame(corMatBrain %>% dplyr::filter(p.value.Control<0.05 & estimate.Control>0) 
                            %>% dplyr::select(gene) %>% mutate(Type="Control"))))

#Negatively correlated between RRBS and RNA-Seq
negCorBrain<-unique(rbind(as.data.frame(corMatBrain %>% dplyr::filter(p.value.IR<0.05 & estimate.IR<0) 
                            %>% dplyr::select(gene) %>% mutate(Type="IR")),
              as.data.frame(corMatBrain %>% dplyr::filter(p.value.Control<0.05 & estimate.Control<0) 
                            %>% dplyr::select(gene) %>% mutate(Type="Control"))))
options(repr.plot.width=6, repr.plot.height=8)
plot_grid(createUpsetPlot(posCorBrain), createUpsetPlot(negCorBrain),ncol=1)
posCorBrain$tissue<-"BRN"
negCorBrain$tissue<-"BRN"

Using Type as value column: use value.var to override.

Using Type as value column: use value.var to override.
In [195]:
allPosCor<-rbind(negCorBrain,negCorRetina)
allPosCor$Type<-paste(allPosCor$tissue,allPosCor$Type,sep="_")

    df<-reshape2::dcast(allPosCor, gene~Type, length)
    Groups = colnames(df)[2:5]
    plotUpset2<-upset(df, Groups, name='Groups', width_ratio=0.1, min_size=0, min_degree=1, set_sizes=FALSE,
      queries=list(
        upset_query(set='BRN_Control', fill='lightblue'),
        upset_query(set='BRN_IR', fill='orange'),
        upset_query(set='RTN_Control', fill='blue'),
        upset_query(set='RTN_IR', fill='brown')
      ),
      themes=upset_default_themes(text=element_text(size=24,face="bold"), 
                                  axis.title=element_text(face="bold",size=26),
                                  panel.border = element_rect(colour = "black", fill=NA, size=1)),
      sort_intersections_by='cardinality',
      base_annotations=list('Size'=(intersection_size(counts=TRUE))),
       matrix=intersection_matrix(geom=geom_point(shape='square filled',size=5)))
options(repr.plot.width=15, repr.plot.height=7)
plotUpset2
dfSum = df %>% mutate(sum = select(., !starts_with("gene") & !contains("Control")) %>% rowSums(na.rm = TRUE))

Using tissue as value column: use value.var to override.
A data.frame: 1 × 5
geneBRN_ControlBRN_IRRTN_ControlRTN_IR
<chr><int><int><int><int>
10610040J01Rik0100
In [200]:
subset(dfSum, sum>1 & BRN_Control==0 & RTN_Control==0)
A data.frame: 9 × 6
geneBRN_ControlBRN_IRRTN_ControlRTN_IRsum
<chr><int><int><int><int><dbl>
165033430I15Rik01012
339Eln 01012
394Fgf11 01012
401Fktn 01012
413G3bp2 01012
454Gm28778 01012
668Ncs1 01012
797Ppm1e 01012
924Sgsh 01012
In [198]:
options(repr.plot.width=20, repr.plot.height=8)
plot_grid(plotUpset,plotUpset2,nrow=1,labels=c("A","B"))
In [203]:
setwd("/home/jupyter/glds202_203/Pathways")
cases<-c("IR","Control")
for (j in 1:length(cases)){
    title<-cases[j]
    tmpData<-subset(as.data.frame(corMatRetina),select=c(paste0("estimate.",title),paste0("p.value.",title),"gene"))
    colnames(tmpData)<-c("estimate","pvalue","gene")
    tmpData<-subset(tmpData,!is.na(estimate))
    tmpData$pvalue<-ifelse(tmpData$pvalue==0, 0.0000000001, tmpData$pvalue)
    #tmpData$stat<-sign(as.numeric(tmpData$estimate))* -log10(as.numeric(tmpData$pvalue))
    tmpData$stat<-as.numeric(tmpData$estimate) * -log10(as.numeric(tmpData$pvalue))
    tmpData<-subset(tmpData,select=c("stat","gene"))
    tmpData<-subset(tmpData,!is.na(stat))

    rankMetric <- as.numeric(tmpData$stat)
    names(rankMetric) <- tmpData$gene
    rankMetric <- sort(rankMetric, decreasing = TRUE)
    
    ####GO####
    gseaGO <- gseGO(geneList=rankMetric,
                                     ont ="BP",
                                     keyType="SYMBOL",
                                     pvalueCutoff = 0.25,
                                     by = "fgsea",
                                     pAdjustMethod = "BH",
                                     OrgDb = org.Mm.eg.db)
    if(dim(gseaGO)[1]>0){
            gseaRes<-data.frame(gseaGO)
            gseaRes$Dataset<-"RTN"
            gseaRes$Group<-title
            write.table(gseaRes,paste0("RTN",title,"_clusterProfiler.GO.BP.all.symbol.tsv"),sep="\t",quote=FALSE,row.names=F)
            }
    }

preparing geneSet collections...

GSEA analysis...

Warning message in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, :
“There are duplicate gene names, fgsea may produce unexpected results.”
no term enriched under specific pvalueCutoff...

preparing geneSet collections...

GSEA analysis...

Warning message in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, :
“There are duplicate gene names, fgsea may produce unexpected results.”
no term enriched under specific pvalueCutoff...
In [205]:
corTissue<-rbind(posCorRetina %>% mutate(Sign="pos"), negCorRetina %>% mutate(Sign="neg"),
                 posCorBrain %>% mutate(Sign="pos"), negCorBrain %>% mutate(Sign="neg"))

corAssay<-rbind(posCorMeth %>% mutate(Sign="pos"), negCorMeth %>% mutate(Sign="neg"),
                 posCorRna %>% mutate(Sign="pos"), negCorRna %>% mutate(Sign="neg"))
corTissue[1,]
subset(ddply(unique(subset(corTissue,Type=="IR",select=c("gene","tissue"))), .(gene), nrow),V1>1)$gene
subset(ddply(unique(subset(corTissue,Type=="Control",select=c("gene","tissue"))), .(gene), nrow),V1>1)$gene
A data.frame: 1 × 4
geneTypetissueSign
<chr><chr><chr><chr>
11700029J07RikIRRTNpos
  1. '5033430I15Rik'
  2. 'Arrdc3'
  3. 'Atad2b'
  4. 'Avpi1'
  5. 'Eln'
  6. 'Fbxo21'
  7. 'Fgf11'
  8. 'Fktn'
  9. 'G3bp2'
  10. 'Gabrg3'
  11. 'Gm28778'
  12. 'Il20ra'
  13. 'Lrrc28'
  14. 'Mfsd9'
  15. 'Mgmt'
  16. 'Mid1'
  17. 'Mrc2'
  18. 'Napa'
  19. 'Ncs1'
  20. 'Nsmce1'
  21. 'Ppm1e'
  22. 'Rmc1'
  23. 'Rpl8'
  24. 'Sgsh'
  25. 'Stim1'
  26. 'Sufu'
  27. 'Tmem39b'
  28. 'Tufm'
  29. 'Uba3'
  30. 'Ulk3'
  31. 'Wrap53'
  32. 'Yipf5'
  33. 'Zfp703'
  34. 'Zfp787'
  1. '1700030C10Rik'
  2. 'Abracl'
  3. 'Adam23'
  4. 'Adar'
  5. 'Agbl2'
  6. 'Anapc15'
  7. 'Btc'
  8. 'Ctnnal1'
  9. 'Exosc8'
  10. 'Fbxo6'
  11. 'Gatd1'
  12. 'Gm4673'
  13. 'Gnat1'
  14. 'Hdgf'
  15. 'Itga8'
  16. 'Kcnip1'
  17. 'Kcnk1'
  18. 'Lrp1b'
  19. 'Mak16'
  20. 'Mta3'
  21. 'Pex11a'
  22. 'Phactr1'
  23. 'Phc2'
  24. 'Pik3cd'
  25. 'Ppcs'
  26. 'Sipa1l3'
  27. 'Socs2'
  28. 'St6galnac2'
  29. 'Tlk1'
  30. 'Tmem18'
  31. 'Trappc6b'
  32. 'Zbtb46'
  33. 'Zfp768'
In [206]:
diffVarBrain<-fread("/home/jupyter/glds202_203/DiffVar/BRN_IR_vs_Control.tsv",header=T,sep="\t")
sort(subset(diffVarBrain,
       external_gene_name %in% c(subset(negCorBrain,Type=="IR")$gene,subset(posCorBrain,Type=="IR")$gene) & 
       (FDR.dispersion<0.05 | FDR.mean<0.05))$external_gene_name)
  1. 'Cox5b'
  2. 'Lamc1'
In [ ]:
tmp1<-subset(counts202,grepl("Cul1",external_gene_name))
tmp2<-subset(counts203,grepl("Cul1",external_gene_name))
In [ ]:
options(repr.plot.width=10, repr.plot.height=8)
tmp1Long<-melt(subset(tmp1,select=-c(ensembl_gene_id))) %>%
            mutate(group=ifelse(grepl("HLLC_IRC",variable),"Control",ifelse(grepl("HLLC_IR_",variable),"IR",
                                                                            ifelse(grepl("HLU_IRC",variable),"HLU","Combined"))),
                  timepoint=ifelse(grepl("4mon",variable),"4m",ifelse(grepl("1mon",variable),"1m","7d")))
ggplot(tmp1Long) + geom_boxplot(aes(x=group,y=value,color=group)) + xlab("gene") + ylab("Normalized counts") + 
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + facet_wrap(.~external_gene_name,scale="free") + ggtitle("Brain")

tmp2Long<-melt(subset(tmp2,select=-c(ensembl_gene_id))) %>%
            mutate(group=ifelse(grepl("HLLC_IRC",variable),"Control",ifelse(grepl("HLLC_IR_",variable),"IR",
                                                                            ifelse(grepl("HLU_IRC",variable),"HLU","Combined"))),
                  timepoint=ifelse(grepl("4mon",variable),"4m",ifelse(grepl("1mon",variable),"1m","7d")))
ggplot(tmp2Long) + geom_boxplot(aes(x=group,y=value,color=group)) + xlab("gene") + ylab("Normalized counts") + 
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + facet_wrap(.~external_gene_name,scale="free") + ggtitle("Retina")

Warning message in melt(subset(tmp1, select = -c(ensembl_gene_id))):
“The melt generic in data.table has been passed a data.frame and will attempt to redirect to the relevant reshape2 method; please note that reshape2 is deprecated, and this redirection is now deprecated as well. To continue using melt methods from reshape2 while both libraries are attached, e.g. melt.list, you can prepend the namespace like reshape2::melt(subset(tmp1, select = -c(ensembl_gene_id))). In the next version, this warning will become an error.”
Using external_gene_name as id variables

Warning message in melt(subset(tmp2, select = -c(ensembl_gene_id))):
“The melt generic in data.table has been passed a data.frame and will attempt to redirect to the relevant reshape2 method; please note that reshape2 is deprecated, and this redirection is now deprecated as well. To continue using melt methods from reshape2 while both libraries are attached, e.g. melt.list, you can prepend the namespace like reshape2::melt(subset(tmp2, select = -c(ensembl_gene_id))). In the next version, this warning will become an error.”
Using external_gene_name as id variables
In [259]:
diffVarRetina<-fread("/home/jupyter/glds202_203/DiffVar/RTN_IR_vs_Control.tsv",header=T,sep="\t")
diffVarRetina<-subset(diffVarRetina,select=c("external_gene_name"), FDR.dispersion<0.1)

for (i in 1:dim(subset(manifest))[1]){
    tmp<-subset(merged203,select=c(subset(manifest)$Retina_rna[i],
                                   paste0("ratio_",subset(manifest)$Retina_meth[i])),
                     external_gene_name %in% diffVarRetina$external_gene_name)
    colnames(tmp)<-c("Brain","Retina")
    if(i==1){
        longRna<-tmp
    }
    else{
        longRna<-rbind(longRna,tmp)
    }
}
ggplot(longRna,aes(x=Brain,y=Retina)) + geom_point(size=5) + geom_smooth(method="lm")

`geom_smooth()` using formula 'y ~ x'
In [233]:
manifest[1,]
A data.table: 1 × 13
SubjectRepTissueSampleTimepointGroupBrain_methRetina_methBrain_rnaRetina_rnaCountBrain_ageRetina_age
<chr><chr><chr><chr><chr><chr><chr><chr><chr><chr><int><dbl><dbl>
M83Rep1BothHLLC_IRC_4mon_Rep1_M834mControlCFG1778CFG1837Mmus_C57_6J_BRN_HLLC_IRC_4mon_Rep1_M83Mmus_C57_6J_RTN_HLLC_IRC_4mon_Rep1_M83413.33.4
In [224]:
#Network of genes with correlation and differential expression/variability
library("STRINGdb")
string_db <- STRINGdb$new(version="11.5", species=10090,score_threshold=400,input_directory="")

diffVarBrain<-fread("/home/jupyter/glds202_203/DiffVar/BRN_IR_vs_Control.tsv",header=T,sep="\t")
diffVarBrain<-subset(diffVarBrain,select=c("external_gene_name"), FDR.dispersion<1 | FDR.mean<1)
diffVarBrain<-unique(subset(diffVarBrain,
       external_gene_name %in% c(subset(posCorBrain,Type=="Control")$gene,
                                subset(negCorBrain,Type=="Control")$gene,
                                subset(posCorBrain,Type=="IR")$gene,
                                subset(negCorBrain,Type=="IR")$gene),
    select=c("external_gene_name")))
length(diffVarBrain$external_gene_name)
cat(paste0(sort(diffVarBrain$external_gene_name),sep=","))
cat("\n")

diffVarRetina<-fread("/home/jupyter/glds202_203/DiffVar/RTN_IR_vs_Control.tsv",header=T,sep="\t")
diffVarRetina<-subset(diffVarRetina,select=c("external_gene_name"), FDR.dispersion<1 | FDR.mean<1)
diffVarRetina<-unique(subset(diffVarRetina,
       external_gene_name %in% c(subset(posCorRetina,Type=="Control")$gene,
                                subset(negCorRetina,Type=="Control")$gene,
                                subset(posCorRetina,Type=="IR")$gene,
                                subset(negCorRetina,Type=="IR")$gene),
    select=c("external_gene_name")))
length(diffVarRetina$external_gene_name)
cat(paste0(sort(diffVarRetina$external_gene_name),sep=","))

df<-merge(diffVarBrain,diffVarRetina,by=c("external_gene_name"))
colnames(df)<-c("gene")

df<-diffVarBrain
colnames(df)<-c("gene")
diffVarBrain$group<-"BRN_IR"
diffVarRetina$group<-"RTN_IR"
df<-rbind(diffVarBrain,diffVarRetina)
1127 
0610040J01Rik, 1500009L16Rik, 1700030C10Rik, 1700034H15Rik, 1700102P08Rik, 2410018L13Rik, 2510009E07Rik, 2700049A03Rik, 2810455O05Rik, 2810459M11Rik, 3010001F23Rik, 3110009E18Rik, 4632411P08Rik, 4930447C04Rik, 4933427D14Rik, 5033430I15Rik, 5330434G04Rik, 5830408C22Rik, 6430590A07Rik, 9030624G23Rik, 9130024F11Rik, 9330111N05Rik, 9330188P03Rik, 9530068E07Rik, A330102I10Rik, A630072M18Rik, A930002C04Rik, Aagab, Abca5, Abcb1b, Abcb4, Abcb8, Abcf1, Abcg1, Abracl, Acadvl, Ache, Ackr3, Acot1, Acot13, Acot4, Acpp, Actn4, Actr1a, Actr2, Adam11, Adam23, Adamtsl5, Adar, Adgra2, Ado, Adra1d, Adra2c, Adssl1, Afg1l, Afg3l2, Afmid, Agbl2, Ago4, AI837181, Akap13, Akr7a5, Alcam, Alg3, Alkal2, Alkbh1, Alkbh5, Alox12b, Anapc15, Anapc2, Ankrd13a, Ankrd13b, Ankrd40, Ankrd45, Ap1s2, Ap3m2, Ap4e1, Ap5s1, Apbb1, Apcdd1, Aplp1, Arcn1, Arhgap28, Arhgap42, Arhgef10l, Arhgef33, Arhgef4, Arhgef5, Arih2, Arl10, Arl6ip4, Arl6ip5, Arl8b, Arrdc3, Arxes1, Arxes2, Asb8, Asnsd1, Ass1, Atad2b, Atg16l2, Atp2b2, Atp5g2, Atp5o, Atpaf1, AU040320, Avpi1, B230219D22Rik, B4galnt4, B4gat1, B630019K06Rik, Bahcc1, Baz2a, Bbs4, Bbs5, Bbs7, BC029722, BC037039, Bcap29, Bcat2, Bcdin3d, Bcl7c, Bean1, Bend4, Bend5, Bicdl1, Bin2, Birc6, Blvra, Bmpr1b, Boc, Brd1, Brd7, Bri3, Brms1l, Bsn, Btbd8, Btc, C030037D09Rik, C1ql3, C2, C230014O12Rik, C2cd2l, C330011M18Rik, Cactin, Cadps, Camk2d, Camk2n2, Caml, Camsap2, Capns1, Capza2, Car11, Car4, Caskin1, Cbr2, Ccdc113, Ccdc158, Ccdc160, Ccdc162, Ccdc22, Ccdc88c, Ccdc9, Ccdc9b, Ccl27a, Ccn1, Ccnk, Cdc14a, Cdc16, Cdc40, Cdh4, Cdh9, Cdk13, Cdk20, Cdk5rap2, Cdkn1b, Cdkn1c, Cdkn2d, Cds2, Cdv3, Celsr3, Cep63, Cep70, Cerox1, Cert1, Cfap298, Cfap36, Cfap418, Cfap57, Cfap61, Cfap97d2, Cftr, Chchd3, Chka, Chordc1, Chst10, Chst15, Chsy3, Chtop, Cilk1, Cipc, Ckap2l, Clcn4, Clcn5, Cldn1, Cldn12, Clspn, Cmtr2, Cnih1, Cnst, Cntnap2, Cntrl, Col9a2, Coq10b, Coq7, Cox17, Cox5b, Cox6c, Cpeb1, Cpne5, Cpq, Cpsf7, Crb2, Crem, Crtc1, Csf1r, Csnk1g3, Ctf1, Ctnnal1, Ctnnb1, Ctsh, Ctsz, Cul5, Cutc, Cxcl12, Cyb5a, Cycs, Cygb, Cyp46a1, Cyth3, D130020L05Rik, D430036J16Rik, Daam2, Dach1, Dcaf10, Dcakd, Dcbld1, Ddit4l, Ddx46, Ddx49, Dennd4c, Desi1, Dhtkd1, Dhx29, Diablo, Dido1, Dip2b, Dip2c, Dis3, Dlat, Dlg4, Dlg5, Dmpk, Dnaaf3, Dnah9, Dnaja3, Dnajb2, Dnajc15, Dnajc3, Dnmbp, Dnmt3a, Dok7, Dolpp1, Dpp6, Dpy19l1, Dpysl3, Dpysl5, Drd1, Dtx4, Dusp16, Dusp27, Dusp8, Dvl3, Dym, Dyrk1b, Ecpas, Edem2, Eef1g, Efcab12, Eid2b, Eif1b, Eif3e, Eif3f, Eif3l, Eif5a, Eif5a2, Eln, Emc10, Eml1, Emp3, En2, Epm2a, Epm2aip1, Eps8, Erbb3, Erich5, Errfi1, Esyt1, Ets1, Etv1, Exoc3l, Exoc3l4, Exosc3, Exosc5, Exosc8, Fads2, Faf2, Fahd1, Faim2, Fam107a, Fam110a, Fam117b, Fam118b, Fam131c, Fam136a, Fam163b, Fam171a1, Fam172a, Fam189a1, Fam207a, Fam234b, Fam53c, Fbln1, Fbn2, Fbxl17, Fbxl20, Fbxl3, Fbxo21, Fbxo6, Fbxw7, Fbxw9, Fdps, Fem1b, Fgd3, Fgd4, Fgf11, Fgf9, Fgfbp3, Fhl3, Fign, Fkbp1b, Fktn, Flcn, Fn3krp, Frk, Fsd1, Fyb, Fzd7, G3bp2, Gabrg3, Gadd45gip1, Galc, Galnt7, Galnt9, Garnl3, Gas7, Gas8, Gatad2b, Gatd1, Gatm, Gcdh, Get3, Gfer, Ggn, Ggta1, Ghsr, Gimap1, Gins4, Glra1, Glt1d1, Gltp, Gm10545, Gm11613, Gm11642, Gm12216, Gm12743, Gm13889, Gm14325, Gm15609, Gm17435, Gm20387, Gm26559, Gm26871, Gm27003, Gm28050, Gm28778, Gm29394, Gm34552, Gm35394, Gm41760, Gm4285, Gm43254, Gm45605, Gm45823, Gm4673, Gm48623, Gm50394, Gm5602, Gm5784, Gm9828, Gm9856, Gnat1, Gng2, Golgb1, Gorasp1, Gpn1, Gpr137, Gpr139, Gpr3, Gpx7, Grb7, Gria2, Grik3, Grk2, Gsn, Gtf2e2, Gtf2h4, Gyg, Gzmm, H1f5, H3c4, H3f3b, H6pd, Haus1, Haus4, Haus5, Haus7, Hcrtr1, Hdgf, Hdgfl3, Herpud2, Hhip, Hic2, Hid1, Hif1an, Higd1a, Hint2, Hmces, Hmgcr, Hnrnph3, Homez, Hook3, Hsd17b7, Hspa8, Hspb8, Ids, Idua, Igbp1, Igfbp7, Igsf3, Il10ra, Il20ra, Il33, Impact, Ino80, Inpp5a, Ints11, Iqgap1, Ireb2, Isyna1, Itga8, Itga9, Itgb4, Jazf1, Jkamp, Jtb, Kat6b, Katna1, Katnb1, Kbtbd3, Kcnh5, Kcnip1, Kcnk1, Kcnq4, Kcp, Kctd6, Kdelr2, Kdm2a, Kdm8, Keap1, Khk, Kif1b, Kif1c, Kif3a, Kifbp, Kifc3, Kirrel2, Klhl1, Klhl14, Klhl35, Kpnb1, Krt12, Kyat3, L2hgdh, Lamc1, Larp1, Lars2, Lcmt2, Lcorl, Lgals8, Lgr6, Lima1, Limd1, Lin7c, Llgl1, Lman1, Lmtk3, Lonrf3, Loxl3, Lpcat3, Lpin2, Lrfn4, Lrp1b, Lrrc28, Lrrc49, Lrrfip2, Lrrn2, Lsm10, Lzts1, Lzts2, Macir, Mak16, Malat1, Mall, Man2a1, Map2k3, Map4k2, Map7d2, Map9, Mapk3, Mapk8, Mapre2, Marchf3, Marf1, Mast1, Mcc, Mcf2l, Mcfd2, Mcrip2, Mdm1, Med10, Med15, Med23, Meg3, Mettl5, Mfsd14b, Mfsd6, Mfsd9, Mgat5, Mgmt, Mgrn1, Mical2, Mid1, Mier1, Mindy2, Mios, Mipol1, Mkrn1, Mllt6, Mms22l, Morc2a, Mospd1, Moxd1, Mpnd, Mpp2, Mpst, Mrc2, Mroh5, Mrpl3, Msl3, Msra, Mta3, Mtg2, Mtss2, Mtus2, Mus81, Mxd4, Myadml2, Mybl2, Mycl, Myh9, Myo1d, Myo1g, Myo9b, Naa10, Naa50, Naa60, Nalcn, Nanos1, Nanp, Napa, Nav1, Ncan, Ncapd2, Nck1, Ncoa2, Ncoa6, Ncor1, Ncs1, Ndst3, Ndufa1, Ndufa12, Ndufa8, Ndufb5, Ndufs6, Necap2, Neil1, Neurod2, Nexn, Nfat5, Nfatc3, Nfatc4, Nfkbie, Nfya, Nfyb, Nhlrc2, Niban2, Nim1k, Nipsnap2, Noct, Nog, Nopchap1, Nova1, Npas4, Npl, Nploc4, Nptn, Nr3c2, Nr4a2, Nrg4, Nrgn, Nrxn2, Nsd1, Nsmce1, Ntmt1, Nubp1, Nudt16l1, Nup188, Nxn, Oaf, Obscn, Olfm1, Olfml2b, Olfml3, Oprm1, Optn, Osbpl1a, Osbpl3, Osbpl6, Ostf1, Otud1, Ovgp1, Oxsr1, P4ha3, P4hb, Pabpc1, Pacs2, Pafah1b1, Pak5, Pals2, Pamr1, Pank1, Papola, Paqr8, Paqr9, Parp1, Parp14, Parp16, Parp3, Pbxip1, Pcdhga1, Pcdhga9, Pcdhgb4, Pcdhgb6, Pced1a, Pck2, Pcnx2, Pcsk6, Pdap1, Pdcd2, Pde6d, Pde7a, Pdgfb, Pdlim2, Pdp1, Pea15a, Per3, Pex11a, Pex12, Pgam2, Pgpep1, Phactr1, Phactr2, Phc1, Phc2, Phc3, Phf13, Phf23, Phpt1, Pias2, Picalm, Pigq, Pigv, Pik3cd, Pim2, Pinx1, Pip4k2b, Pitpnm1, Pkp2, Pla2g12a, Plce1, Plekha4, Plekhd1, Plekhg5, Plekhh1, Plekhm3, Plppr1, Pls1, Pm20d2, Pnma8b, Pno1, Poglut1, Poldip2, Pole, Poll, Polr3a, Polr3b, Ppa2, Ppara, Ppcs, Pphln1, Ppia, Ppic, Ppip5k1, Ppm1e, Ppm1k, Ppp1r10, Ppp1r16b, Ppp1r3f, Ppp5c, Ppt1, Prcc, Prkab2, Prkcg, Prkch, Prkdc, Prkn, Prps2, Psmb10, Psph, Psrc1, Ptbp1, Ptcd2, Pter, Ptpn1, Ptpn13, Ptpn18, Ptprm, Ptprn, Ptpro, Puf60, Pwp1, Pxdc1, Pxdn, Pxmp2, Qrich1, Qrsl1, Rab27a, Rab33a, Rab40b, Rab42, Rabggta, Rabl6, Rad9a, Radil, Ramp1, Ranbp17, Ranbp3l, Rap1gds1, Rasa4, Rasgrp1, Rassf2, Rassf8, Rbm26, Rbm4, Rcan3, Rcbtb2, Reck, Reep4, Retreg1, Rfxap, Rgcc, Rgs7, Rgs7bp, Rlf, Rmc1, Rmi1, Rnf114, Rnf144a, Rnf166, Rnf170, Rnf182, Rnf185, Rnf217, Rnft2, Robo2, Rogdi, Rpgrip1, Rpia, Rpl12, Rpl15, Rpl24, Rpl26, Rpl29, Rpl34, Rpl35a, Rpl41, Rpl7, Rpl8, Rplp0, Rplp2, Rpn1, Rps6ka3, Rrbp1, Rreb1, Rrp8, Rrs1, Rsc1a1, Rwdd3, Rxfp3, S100a11, S1pr5, Sae1, Samd11, Samd15, Sart1, Scfd2, Scly, Scn2b, Scnn1a, Sdf2, Sema6a, Sema6c, Septin10, Serf2, Sergef, Serhl, Serpine2, Serpinh1, Sesn3, Setd4, Sgcd, Sgsh, Sh3bgrl2, Sh3glb1, Shank1, Shc3, Shc4, Shd, Shf, Shisa6, Shisa8, Siae, Sik2, Sipa1l3, Slc17a6, Slc1a4, Slc24a2, Slc24a3, Slc25a15, Slc25a20, Slc25a5, Slc2a8, Slc35b2, Slc35e2, Slc38a2, Slc38a6, Slc39a9, Slc41a2, Slc41a3, Slc44a1, Slc47a1, Slc66a2, Slc6a11, Slc7a1, Slc8a3, Slc9a1, Slco2a1, Slco2b1, Slco4a1, Slk, Slu7, Smad1, Smarcc2, Smc1a, Smc4, Smc6, Smg7, Smg8, Smurf1, Smurf2, Snhg4, Snph, Snrpd3, Snx32, Socs2, Socs3, Soga1, Sox7, Spaar, Spag4, Specc1l, Src, Srcap, Srek1, Srek1ip1, Srp54a, Srp72, Srrm4, Ssb, Ssbp3, Sst, Ssx2ip, St6gal2, St6galnac2, Stim1, Stim2, Stip1, Stk32b, Stpg1, Strbp, Strn4, Stxbp4, Stxbp5l, Sucla2, Sufu, Sugp2, Sult4a1, Surf2, Swt1, Sympk, Syngr2, Synpo2l, Syt5, Tafa2, Tamm41, Tank, Tbc1d2, Tbc1d22a, Tbc1d30, Tc2n, Tceanc2, Tcerg1, Tecpr1, Tenm2, Terb1, Tesk2, Tet3, Tfam, Tgm2, Thap11, Thoc5, Thoc6, Thop1, Tigar, Tjap1, Tjp1, Tjp2, Tlcd2, Tlcd3a, Tlk1, Tm2d1, Tmc3, Tmcc2, Tmco4, Tmed5, Tmem108, Tmem160, Tmem18, Tmem204, Tmem250-ps, Tmem35a, Tmem35b, Tmem37, Tmem38b, Tmem39a, Tmem39b, Tmem59l, Tmie, Tmpo, Tmtc4, Tnfsfm13, Tnks, Togaram1, Tom1l1, Tomm40l, Top2a, Tor1aip2, Tox2, Tpbgl, Tpd52l1, Tppp3, Trafd1, Trappc2l, Trappc6b, Trim62, Trim71, Trim9, Trip11, Trip12, Trmt2a, Trnp1, Trpc5, Trpm3, Ttc19, Ttc3, Ttpa, Tubd1, Tufm, Tyw3, Uba3, Ubald2, Ube2b, Ube2cbp, Ube2g1, Ube2n, Ube3c, Ubxn2b, Uchl3, Ucp2, Ufc1, Uhrf1bp1l, Ulk3, Umps, Unc119, Unc5d, Urm1, Usp15, Usp3, Usp45, Usp5, Ust, Utp15, Utp25, Uvrag, Vamp8, Vasp, Vcpkmt, Vdac1, Vgll3, Vps13a, Vps13b, Vps26a, Vps4a, Vsig10, Vwa8, Wdhd1, Wdr33, Wdr4, Wdr74, Wdr75, Wdr81, Wfs1, Wipf2, Wnk3, Wrap53, Wrap73, Wtip, Xkr4, Xkr6, Ybx2, Yeats2, Yipf5, Ythdf2, Ywhae, Zbtb1, Zbtb10, Zbtb4, Zbtb46, Zcchc18, Zcchc3, Zfat, Zfhx3, Zfp1, Zfp185, Zfp2, Zfp236, Zfp28, Zfp282, Zfp365, Zfp382, Zfp385b, Zfp407, Zfp414, Zfp472, Zfp516, Zfp523, Zfp609, Zfp688, Zfp703, Zfp768, Zfp770, Zfp777, Zfp786, Zfp787, Zfp853, Zfp867, Zfp94, Zfp971, Zfp994, Zfp995, Zfyve26, Zic3, Zmpste24, Znrd2, Zrsr2, Zswim4,
1062 
1110032F04Rik, 1110038B12Rik, 1700029J07Rik, 1700030C10Rik, 1810024B03Rik, 1810059H22Rik, 2610044O15Rik8, 2610306M01Rik, 3110056K07Rik, 4930515G01Rik, 4932441J04Rik, 4933404O12Rik, 4933407K13Rik, 5033430I15Rik, 5330439K02Rik, 6430571L13Rik, 6720427I07Rik, 9930014A18Rik, A230072C01Rik, A230077H06Rik, A430035B10Rik, A830052D11Rik, A930012O16Rik, Abat, Abca13, Abhd5, Ablim3, Abracl, Abt1, Abtb1, Acbd6, Ache, Ackr3, Actr1a, Acvrl1, Adam23, Adamts16, Adamts19, Adamtsl3, Adar, Adarb1, Adck2, Adck5, Adprm, Adra2b, Agap1, Agbl2, Ago3, Agpat2, Aif1, Aig1, Ak4, Akap8l, Aldh18a1, Aldoa, Alkbh8, Aloxe3, Als2cl, Amd1, Anapc1, Anapc15, Angptl4, Ankrd13c, Ankrd33b, Aopep, Ap1s3, Ap3d1, Apba3, Aplp1, App, Araf, Arap2, Arf1, Arf4, Arf5, Arfgef1, Arg2, Arhgap12, Arhgap19, Arhgap30, Arhgap33, Arhgef15, Arid1b, Arid5a, Arl14ep, Arpp21, Arrdc3, Ascc2, Ascl1, Ash2l, Asic1, Asl, Asphd1, Asxl1, Atad2b, Atat1, Atg2a, Atp2b4, Atp5a1, Atp5g3, Atp5pb, Atp6ap1, Atp6v1a, Atp6v1c1, Atp6v1h, Atp9a, Avpi1, B230206L02Rik, B4galt7, B9d1, Babam1, Bahcc1, Baz1b, BC002059, BC005537, BC034090, BC037704, Bcas3, Bcl3, Bdp1, Bend5, Bfar, Bfsp2, Birc6, Blmh, Bloc1s3, Bmp6, Bmt2, Braf, Brca2, Brd3os, Brinp1, Btbd10, Btbd3, Btbd9, Btc, C2, Cacng8, Calcr, Calhm2, Caly, Camk2a, Cant1, Capn1, Card19, Carm1, Casz1, Cav1, Cbr3, Cbx2, Ccdc116, Ccdc117, Ccdc137, Ccdc84, Cchcr1, Cckbr, Ccnd1, Cd276, Cd99l2, Cdc14a, Cdc25a, Cdc25b, Cdk2ap1, Cdk5r1, Cdk5r2, Cdkn1a, Cdon, Cebpd, Cebpzos, Cela1, Celf3, Celf4, Cep135, Cfap36, Cfap53, Cgrrf1, Chchd6, Chdh, Chil1, Chmp2a, Chmp3, Chmp5, Chmp6, Chp1, Chrm3, Chrnb3, Chst14, Chst2, Chsy1, Ciao3, Cic, Cipc, Clcn3, Cldn3, Clip1, Clk2, Clstn1, Cnksr3, Coasy, Cobl, Colec12, Cops3, Cops7b, Cops9, Coq8a, Coro1b, Cox6c, Cplx1, Cracdl, Crb1, Crispld2, Crk, Cry2, Csnk1d, Csnk1g2, Csnk2a1, Cspg4, Ctbp1, Ctnnal1, Ctsa, Cttnbp2, Ctu2, Ctxn1, Cyb5d2, Cyth4, D16Ertd472e, D430040D24Rik, D630045J12Rik, D830050J10Rik, Daam1, Dbndd1, Dctpp1, Dcun1d5, Dcx, Ddi2, Ddn, Ddost, Ddr2, Ddx3x, Ddx5, Ddx6, Degs1, Degs2, Dennd4b, Derl3, Dgcr8, Dgkq, Dgkz, Dhx34, Dhx57, Disp1, Dlec1, Dleu2, Dnaaf9, Dnajc12, Dnajc25, Dnajc6, Dnttip1, Dock8, Dock9, Dot1l, Dpf1, Dusp16, Dusp3, Dync1i2, Dync2i2, E130114P18Rik, Ebf2, Eci1, Ecscr, Edn3, Efcab15, Efhd2, Efs, Egln3, Eid2b, Eif2ak4, Eif2b1, Eif2d, Elf4, Elmo3, Eln, Elovl2, Elovl7, Emc9, Emg1, Endou, Engase, Enoph1, Enpp1, Entpd7, Eogt, Epb41l4b, Epha2, Epha6, Epor, Eps8l1, Etnk1, Evc, Evi5, Exosc8, Eya4, Ezh2, F630040K05Rik, Fam102a, Fam104a, Fam149b, Fam193b, Fam20b, Fam210b, Fam214a, Fam234a, Fam89a, Fasn, Fbl, Fbln7, Fbxl14, Fbxl17, Fbxo21, Fbxo31, Fbxo34, Fbxo41, Fbxo6, Fcgr2b, Fdx1, Fes, Fgf11, Fh1, Ficd, Fkbp11, Fkbp14, Fkbp15, Fktn, Flywch2, Fmnl1, Fmnl3, Fmod, Fmr1, Fndc4, Foxn2, Foxn3, Foxp1, Frat1, Fstl4, Ftsj1, Fut10, Fut8, Fv1, Fxyd6, G3bp2, Gab1, Gabra1, Gabrg3, Gal3st3, Galc, Galnt1, Galnt10, Galnt15, Gapvd1, Garem1, Gas2l3, Gatd1, Gdap1l1, Gdf10, Gfy, Ghitm, Gins2, Gipc1, Gkap1, Glrb, Glrx3, Gm11175, Gm12905, Gm13205, Gm14399, Gm14964, Gm15672, Gm19345, Gm19410, Gm20109, Gm20421, Gm26575, Gm26592, Gm26777, Gm28778, Gm33195, Gm38431, Gm40841, Gm42742, Gm4673, Gm49066, Gm6712, Gm9903, Gm9925, Gmcl1, Gmnn, Gnat1, Gnpda1, Gorasp2, Gosr1, Gpc5, Gpkow, Gprasp1, Grb7, Grik5, Grin1, Grk4, Gsx2, Gtf2f2, Gtf3c6, H1f1, H2-T23, H3f3a, Haus3, Hdgf, Hexim2, Hhex, Hilpda, Hint1, Hltf, Hmgb3, Hmox1, Hnrnpr, Hpf1, Hps4, Hsd17b12, Hsp90b1, Hspg2, Hyi, Ice2, Idh3a, Ifitm1, Ift140, Igf2r, Igfbp4, Igfbp6, Ikzf1, Il10rb, Il15ra, Il17d, Il17rc, Il20ra, Inca1, Ing1, Ing4, Inpp5j, Ints11, Irf4, Ist1, Itga2, Itga8, Itpa, Itpka, Itsn2, Jakmip1, Jrk, Jun, Junb, Katnbl1, Kazald1, Kcna5, Kcnc1, Kcng3, Kcnip1, Kcnj11, Kcnj3, Kcnk1, Kcnma1, Kdm6b, Kl, Klc1, Klf3, Klf4, Klf5, Klhl11, Klhl32, Kmt5a, Kyat1, Large1, Ldlrap1, Letmd1, Lfng, Limch1, Lipt2, Lmod1, Lmtk2, Lpcat4, Lrba, Lrif1, Lrp1b, Lrrc28, Lrrc45, Lrrc7, Lrrc75a, Luc7l, Lypd1, Lypd6b, Lypla2, Macrod1, Maea, Magee1, Mak16, Maml1, Map10, Map3k15, Map3k7, Map6, Marchf7, Marchf8, Mat2a, Mbtd1, Mbtps1, Mbtps2, Mchr1, Mcm6, Mecom, Med12l, Med9, Mef2d, Metap2, Mettl1, Mettl27, Mfrp, Mfsd4a, Mfsd9, Mgat4b, Mgmt, Micu2, Mid1, Mief1, Miga2, Mir8114, Mlxipl, Mmd, Mmp17, Mmp28, Mogs, Morc2b, Mrc2, Mrpl24, Mrpl34, Mrpl58, Mrps30, Mrps31, Msh2, Mta1, Mta3, Mtch2, Mtfp1, Mtg1, Mthfr, Mtif3, Mtor, Mtrfr, Myg1, Myl12a, Myo15, Myo5c, Myo9a, N4bp2l1, Naa30, Naa40, Nadk2, Nampt, Napa, Nars, Nat8l, Nbeal2, Ncam1, Ncaph2, Ncdn, Nceh1, Ncs1, Nde1, Ndor1, Ndufa6, Ndufb2, Necab1, Nek7, Nes, Net1, Nfat5, Nfil3, Nfrkb, Nfya, Nhlh2, Nhlrc1, Nhs, Nkrf, Nln, Nnt, Nop2, Nop56, Npas3, Nphs1, Npr3, Nprl3, Nrde2, Nrip1, Nrn1l, Nsmce1, Nub1, Nudcd1, Nudt18, Nudt22, Nup155, Nup210l, Obi1, Ocrl, Odad3, Odad4, Odf2l, Ogfrl1, Olfml2b, Orc4, Oscp1, Oser1, Otud7b, Otx2, P3h3, Pa2g4, Pafah2, Palb2, Palld, Palm3, Panx1, Papss2, Park7, Parva, Pax2, Pax6, Pcbp3, Pcdh15, Pcdha3, Pcdhga2, Pcgf3, Pcmt1, Pcmtd1, Pde1c, Pde6h, Pde8b, Pdf, Pdia6, Pdk2, Pdlim4, Pdpk1, Per1, Pex11a, Pfkfb4, Pfkp, Pgap4, Pgk1, Pgm3, Phactr1, Phc2, Phf5a, Phkg2, Phospho2, Phyh, Pi4ka, Pi4kb, Pigs, Pik3c2b, Pik3cd, Pik3r4, Pir, Pitpna, Pkig, Pknox2, Plekhb2, Plekhf1, Plekhf2, Plekhm3, Plekho1, Plk2, Plpp2, Pltp, Pmpca, Pnpla2, Podxl, Pola1, Pola2, Poldip3, Polq, Polr1a, Polr3d, Polrmt, Porcn, Ppcs, Ppfibp1, Ppm1e, Ppp1cc, Ppp1r13l, Ppp1r1a, Ppp1r3d, Ppp1r9a, Ppp4r3b, Pramel12, Prcp, Prdm16, Prep, Prkaa1, Prkar2b, Prkcsh, Prkd2, Prkd3, Prox1, Psd2, Psen1, Psma4, Psmc1, Psmd1, Ptbp3, Ptdss2, Ptgs1, Ptms, Ptpa, Ptpn13, Ptpn14, Ptprs, Pum3, Purb, Pxmp4, Qdpr, Rab10os, Rab11fip5, Rab12, Rab13, Rab14, Rab33b, Rab36, Rab3il1, Rabep2, Rad52, Raf1, Rap2c, Raph1, Rasgef1c, Rb1cc1, Rbbp7, Rbbp8, Rbfox1, Rbm25, Rbm26, Rbmx, Rbpms, Rcl1, Rfc3, Rgcc, Rgl1, Rhbdd3, Rheb, Rhno1, Rhof, Ric3, Ripk2, Ripor1, Rmc1, Rnf103, Rnf111, Rnf126, Rnf135, Rnf139, Rnf146, Rnf149, Rnf17, Rnf220, Rnh1, Rnpep, Robo4, Rpap3, Rpl10, Rpl19, Rpl21, Rpl32, Rpl7l1, Rpl8, Rprm, Rps11, Rps26, Rps6ka2, Rps6ka3, Rras, Rreb1, Rrm1, Rsbn1l, Rtraf, Rufy2, Rwdd1, Rxrb, S1pr1, S1pr2, Sbsn, Sc5d, Scamp4, Scara3, Scn4b, Scrn1, Scyl3, Sdf4, Sdhd, Sec62, Sema4g, Senp3, Septin1, Serf1, Serinc2, Serpinb6a, Serping1, Setd1b, Sf3b2, Sf3b4, Sgsh, Sgsm1, Sgsm3, Sh2b1, Sh3bgrl, Sh3bp4, Sh3kbp1, Sh3pxd2b, Sh3rf1, Shc2, Shc3, Shd, Shisa2, Sik2, Simc1, Sin3a, Sin3b, Sipa1, Sipa1l1, Sipa1l3, Slbp, Slc25a1, Slc25a22, Slc27a3, Slc27a4, Slc30a1, Slc35d1, Slc35g1, Slc43a3, Slc4a3, Slc4a8, Slc5a3, Slc7a4, Slc9a3r1, Slc9a5, Slco2a1, Slco5a1, Slit2, Smarcad1, Smim24, Smoc1, Sms, Snrpn, Snx2, Soat1, Socs2, Soga3, Sorcs1, Sos2, Sox12, Sox9, Sp4, Spag1, Spag9, Spata2, Spata2l, Spc25, Spdl1, Spg20, Spg7, Spns2, Spopl, Spryd3, Spsb4, Spx, Sqle, Srbd1, Srd5a3, Srrm3, Srsf11, Srsf7, Ssr4, Sstr1, Sstr4, St6galnac2, St8sia3, Stac, Stambp, Stambpl1, Stat5a, Stat6, Stc2, Stim1, Stradb, Stx17, Sufu, Sult2b1, Sybu, Syf2, Syt4, Tada2a, Tafa1, Tapbp, Tarbp2, Tbc1d1, Tbc1d30, Tbc1d32, Tbc1d5, Tbce, Tbl2, Tbp, Tbrg1, Tcerg1l, Tcp1, Tef, Tfcp2, Thap12, Tiam1, Timm44, Timm8a1, Tipin, Tjp1, Tle1, Tle5, Tlk1, Tlk2, Tll1, Tlr2, Tmcc3, Tmed9, Tmem11, Tmem132e, Tmem145, Tmem150c, Tmem164, Tmem170b, Tmem18, Tmem214, Tmem229b, Tmem39b, Tmem64, Tmem94, Tmem97, Tmsb10, Tnfrsf10b, Tnfrsf12a, Tnfrsf13c, Tnfrsf23, Tnks1bp1, Tnnt2, Tns2, Tob1, Tor1aip1, Tra2b, Trappc2, Trappc6b, Trim39, Trim46, Trim8, Trp53, Trpm1, Trpv2, Trpv4, Tsc2, Tsen2, Tsen34, Tshz2, Tspan11, Tsr1, Ttc13, Ttc21a, Tti2, Ttll11, Ttll3, Ttll7, Ttn, Tub, Tubb3, Tubg2, Tubgcp5, Tufm, Txlng, Txn2, Txndc12, Tysnd1, Uba3, Uba6, Ubap2, Ubc, Ube2f, Ubr1, Ubxn8, Ugp2, Ugt8a, Ulk2, Ulk3, Upf2, Upf3a, Uqcrc1, Urgcp, Usp37, Usp42, Usp48, Utp3, Utp6, Vdac3, Vegfc, Vps54, Vps9d1, Vsir, Vstm2b, Vsx2, Wasl, Wdfy1, Wdr4, Wdr47, Wfs1, Wipf2, Wnk2, Wnk4, Wrap53, Xk, Xkr4, Yeats2, Yipf5, Yod1, Zbtb17, Zbtb21, Zbtb46, Zc3h7a, Zer1, Zfp105, Zfp106, Zfp142, Zfp157, Zfp174, Zfp267, Zfp335os, Zfp354a, Zfp369, Zfp384, Zfp445, Zfp469, Zfp53, Zfp563, Zfp598, Zfp599, Zfp606, Zfp637, Zfp661, Zfp703, Zfp750, Zfp768, Zfp78, Zfp784, Zfp787, Zfp871, Zfp9, Znrf2, Zrsr1,
In [ ]:
library(ggvenn)
library(RColorBrewer)
genelist<-subset(merge(tmp, subset(ddply(tmp,.(external_gene_name),nrow),V1>1), by=c("external_gene_name")),external_gene_name!="")
x <- list(
      `RTN_HLU+IR`=subset(genelist, group=="RTN_Combined")$external_gene_name,
      BRN_HLU=subset(genelist, group=="BRN_HLU")$external_gene_name,
      `BRN_HLU+IR`=subset(genelist, group=="BRN_Combined")$external_gene_name,
RTN_HLU=subset(genelist, group=="RTN_HLU")$external_gene_name)
options(repr.plot.width=19, repr.plot.height=20, stroke_alpha=0)
ggvenn(x, show_elements = T, label_sep = ",", fill_color = c("lightblue", "yellow", "lightpink", "grey"),
       text_size=6,stroke_color="white") + ggtitle("C")
diffVarRetina<-fread("/home/jupyter/glds202_203/DiffVar/BRN_Combined_vs_Control.tsv",header=T,sep="\t")
subset(corTissue,grepl("Efna5",gene))
subset(diffVarRetina,grepl("Efna5",external_gene_name))
A data.frame: 0 × 4
geneTypetissueSign
<chr><chr><chr><chr>
A data.table: 1 × 10
ensembl_gene_idControlvsCombined.mean.log2FC.0Statistics.meanPvalue.meanFDR.meanControlvsCombined.dispersion.log2FC.0Statistics.dispersionPvalue.dispersionFDR.dispersionexternal_gene_name
<chr><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><chr>
ENSMUSG00000048915-0.9195698-3.2844610.0010217760.06859507-1.749051-1.1993720.23038311Efna5
In [222]:
tmp<-subset(df, group=="BRN_IR")
tmp$color<-ifelse(tmp$external_gene_name %in% subset(negCorBrain,Type=="IR")$gene, "blue",
                 ifelse(tmp$external_gene_name %in% subset(posCorBrain,Type=="IR")$gene, "red", "gray"))
mapped = string_db$map(tmp, "external_gene_name",removeUnmappedRows = TRUE )
hits <- mapped$STRING_id
In [223]:
# post payload information to the STRING server
options(repr.plot.width=10, repr.plot.height=10)
payload_id <- string_db$post_payload(mapped$STRING_id, colors=mapped$color )
string_db$plot_network(hits, payload_id=payload_id )
In [179]:
clustersList <- string_db$get_clusters(mapped$STRING_id)
options(repr.plot.width=10, repr.plot.height=10)
par(mfrow=c(1,1))
for(i in seq(5:5)){
    tmp<-as.data.frame(clustersList[[i]])
    colnames(tmp)<-c("STRING_id")
    tmp<-merge(tmp,mapped)
    payload_id <- string_db$post_payload(tmp$STRING_id, colors=tmp$color )
    string_db$plot_network(tmp$STRING_id, payload_id=payload_id)
 }
In [159]:
plot_grid(ggplot(manifest, aes(x=Group,y=Retina_age)) + geom_boxplot(position=position_dodge(0.8),width=0.5) + theme(legend.position="none") + geom_jitter(aes(color=Subject),size=5,position=position_dodge(0.3)) + ggtitle("Retina"), 
          ggplot(manifest, aes(x=Group,y=Brain_age)) + geom_boxplot(position=position_dodge(0.8),width=0.5) + theme(legend.position="none") + geom_jitter(aes(color=Subject),size=5,position=position_dodge(0.3)) + ggtitle("Brain"))

Warning message:
“Removed 1 rows containing non-finite values (stat_boxplot).”
Warning message:
“Removed 1 rows containing missing values (geom_point).”
Warning message:
“Removed 3 rows containing non-finite values (stat_boxplot).”
Warning message:
“Removed 3 rows containing missing values (geom_point).”
In [203]:
library(pheatmap)
corGeneList<-subset(corMatRetina, p.value.Combined<0.05 & estimate.Combined<0 & grepl("Slc",gene))$gene

tmp1<-t(apply(subset(merged203, select=c(subset(manifest,Group=="Combined")$Retina_rna), external_gene_name %in% corGeneList),1,scale))
tmp2<-t(apply(subset(merged203, select=paste0("ratio_",c(subset(manifest,Group=="Combined")$Retina_meth)), external_gene_name %in% corGeneList),1,scale))

tmpMerged<-as.data.frame(cbind(tmp1,tmp2))
colnames(tmpMerged)<-c(paste0(subset(manifest,Group=="Combined")$Subject,"_",subset(manifest,Group=="Combined")$Retina_rna),
                      paste0(subset(manifest,Group=="Combined")$Subject,"_",subset(manifest,Group=="Combined")$Retina_meth))

labels<-as.data.frame(colnames(tmpMerged))
colnames(labels)<-c("Sample")
labels$Subject<-c(subset(manifest,Group=="Combined")$Subject,subset(manifest,Group=="Combined")$Subject)
labels$Group<-ifelse(grepl("Mmus",labels$Sample),"RNA","RRBS")
labels <- labels %>% column_to_rownames(var="Sample")

tmpMerged<-tmpMerged[ , order(names(tmpMerged))]
rownames(tmpMerged)<-NULL
rownames(tmpMerged)<-paste0(subset(merged203, external_gene_name %in% corGeneList)$external_gene_name,"\n",
                           subset(merged203, external_gene_name %in% corGeneList)$ensembl_gene_id)

tmpMerged[, "sd"] <- apply(tmpMerged, 1, sd)
tmpMerged<-subset(tmpMerged,sd>0)
tmpMerged<-subset(tmpMerged,select= -c(sd))
dim(tmpMerged)
  1. 5
  2. 8
In [205]:
options(repr.plot.width=6, repr.plot.height=5)
pheatmap(tmpMerged,fontsize = 14,annotation_col=labels,
        show_rownames=TRUE,show_colnames=FALSE,cluster_cols=TRUE,
        clustering_distance_rows = "correlation",clustering_distance_cols = "correlation",
        border_color = NA, main="+ve, Combined")
In [174]:
corGeneList<-subset(corMatMeth, p.value.IR<0.05 & estimate.IR<0)$gene

tmp1<-t(apply(subset(mergedMeth, select=paste0("ratio_",c(subset(manifest,Group=="IR")$Brain_meth)), external_gene_name %in% corGeneList),1,scale))
tmp2<-t(apply(subset(mergedMeth, select=paste0("ratio_",c(subset(manifest,Group=="IR")$Retina_meth)), external_gene_name %in% corGeneList),1,scale))

tmpMerged<-as.data.frame(cbind(tmp1,tmp2))
colnames(tmpMerged)<-c(paste0(subset(manifest,Group=="IR")$Subject,"_",subset(manifest,Group=="IR")$Brain_meth),
                      paste0(subset(manifest,Group=="IR")$Subject,"_",subset(manifest,Group=="IR")$Retina_meth))

labels<-as.data.frame(colnames(tmpMerged))
colnames(labels)<-c("Sample")
labels$Subject<-c(subset(manifest,Group=="IR")$Subject,subset(manifest,Group=="IR")$Subject)
labels$Group<-c(rep("BRN",length(subset(manifest,Group=="IR")$Brain_meth)), rep("RTN",length(subset(manifest,Group=="IR")$Retina_meth)))
labels <- labels %>% column_to_rownames(var="Sample")

tmpMerged<-tmpMerged[ , order(names(tmpMerged))]
rownames(tmpMerged)<-NULL
rownames(tmpMerged)<-paste0(subset(mergedMeth, external_gene_name %in% corGeneList)$external_gene_name,"_",
                           subset(mergedMeth, external_gene_name %in% corGeneList)$ensembl_gene_id)

tmpMerged[, "sd"] <- apply(tmpMerged, 1, sd)
tmpMerged<-subset(tmpMerged,sd>0)
tmpMerged<-subset(tmpMerged,select= -c(sd))
dim(tmpMerged)
  1. 213
  2. 10
In [ ]:
options(repr.plot.width=5, repr.plot.height=5)
pheatmap(tmpMerged,fontsize = 16,annotation_col=labels,
        show_rownames=FALSE,show_colnames=FALSE,cluster_cols=FALSE,
        clustering_distance_rows = "correlation",clustering_distance_cols = "correlation",
        border_color = NA, main="-ve, IR")
In [ ]:
simMatrix <- calculateSimMatrix(as.data.frame(oraGO_Combined)$ID,
                                orgdb="org.Mm.eg.db",
                                ont="BP",
                                method="Rel")
    scores <- setNames(-log10(as.data.frame(oraGO_Combined)$`p.adjust`), as.data.frame(oraGO_Combined)$ID)
    reducedTerms_Combined <- reduceSimMatrix(simMatrix,
                                scores,
                                orgdb="org.Mm.eg.db")

ggplot(reducedTerms_Combined, aes(area = score, subgroup = parentTerm, fill = parentTerm, label = term)) + 
  geom_treemap(size=0.5, col='white') + ggtitle("Combined") +
  geom_treemap_text(col='gray',place="center",reflow=TRUE) +
  geom_treemap_subgroup_border(col='white',size=1) +
  geom_treemap_subgroup_text(col='white',place="center",fontface="bold",size=12,reflow=TRUE,grow=TRUE,min.size=1) + theme(plot.title = element_text(hjust=0.5,size=10)) +
  guides(fill=F)

preparing gene to GO mapping data...

preparing IC data...

Warning message in calculateSimMatrix(as.data.frame(oraGO_Combined)$ID, orgdb = "org.Mm.eg.db", :
“Removed 1 terms that were not found in orgdb for BP”
Warning message:
“`guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.”